#ifndef UHM_DISCRETE_EXTENSION_MESH2PART_HXX_
#define UHM_DISCRETE_EXTENSION_MESH2PART_HXX_

// METIS
extern "C" {
  extern void METIS_PartGraphRecursive(Int_ *, Int_ *, Int_ *, Int_ *,
                                       Int_ *, Int_ *, Int_ *, Int_ *,
                                       Int_ *, Int_ *, Int_ *);
}

namespace uhm {

  template<class M_> Bool_ mesh2part(M_ &mesh, Int_ nparts);

  // partition info is all zero
  template<class M_> Bool_ mesh2part(M_ &mesh, Int_ nparts) {
    
    std::vector< typename M_::Element_* > els;
    mesh.dump(els, 0);

    time_in("uhm:: mesh2part",UHM_TIME_LEVEL_INTERNAL);
    if (els.size() > 1) {

      Int_ n=els.size(), wgtflag=1, numflag=0, options=0, edgecut=0;

      std::vector< Int_ > xadj, adjncy, adjwgt, part;
      part.assign(n, 0);
      elements2graph<M_>(&els, &xadj, &adjncy, &adjwgt, 0);

      //graph_show(std::cout, xadj, adjncy, adjwgt);

      // Partition mesh
      METIS_PartGraphRecursive(&n, &xadj[0], &adjncy[0], NULL, &adjwgt[0],
                               &wgtflag, &numflag,
                               &nparts, &options, &edgecut, &part[0]);

      // Change part
      for (Int_ i=0;i<part.size();++i) 
        els.at(i)->id().set_part(part.at(i));
    }
    time_out(UHM_TIME_LEVEL_INTERNAL);

    return true;
  }
}

#endif
